Past climatic refugia and landscape resistance explain spatial genetic structure in Oriental beech in the South Caucasus

Abstract Predicting species‐level effects of climatic changes requires unraveling the factors affecting the spatial genetic composition. However, disentangling the relative contribution of historical and contemporary drivers is challenging. By applying landscape genetics and species distribution modeling, we investigated processes that shaped the neutral genetic structure of Oriental beech (Fagus orientalis), aiming to assess the potential risks involved due to possible future distribution changes in the species. Using nuclear microsatellites, we analyze 32 natural populations from the Georgia and Azerbaijan (South Caucasus). We found that the species colonization history is the most important driver of the genetic pattern. The detected west–east gradient of genetic differentiation corresponds strictly to the Colchis and Hyrcanian glacial refugia. A significant signal of associations to environmental variables suggests that the distinct genetic composition of the Azerbaijan and Hyrcanian stands might also be structured by the local climate. Oriental beech retains an overall high diversity; however, in the context of projected habitat loss, its genetic resources might be greatly impoverished. The most affected are the Azerbaijan and Hyrcanian populations, for which the detected genetic impoverishment may enhance their vulnerability to environmental change. Given the adaptive potential of range‐edge populations, the loss of these populations may ultimately affect the specie's adaptation, and thus the stability and resilience of forest ecosystems in the Caucasus ecoregion. Our study is the first approximation of the potential risks involved, inducing far‐reaching conclusions about the need of maintaining the genetic resources of Oriental beech for a species' capacity to cope with environmental change.


| INTRODUC TI ON
Trees are known to be playing a substantial role in mitigating the effects of climate change (Anderegg et al., 2020;Bastin et al., 2019).
Yet, their long-term resilience and adaptability is dependent upon genetic diversity, which is currently threatened by climate change itself and by anthropogenic losses of trees on a global scale (Alberto et al., 2013;Hoban et al., 2020;Pauls et al., 2013). Maintaining a high level of genetic diversity and connectivity across the landscape should be a conservation priority, particularly in the world's biodiversity hotspots (Bastin et al., 2019;Fady et al., 2016;Trew & Maclean, 2021).
The Caucasus ecoregion (Figure 1) is one of the biologically richest yet most highly anthropogenically threatened area (Mittermeier et al., 2011;Nikolaishvili & Dvalashvili, 2015;Shatberashvili et al., 2016) and is at high risk of climate change, especially in its eastern part (IPCC, 2022;Nikolaishvili & Dvalashvili, 2015;Shatberashvili et al., 2016). An alarming forest cover loss in the South Caucasus region (i.e., Georgia, Azerbaijan and Armenia) is predicted in this century due to the climatic crisis (Dagtekin et al., 2020;Dering et al., 2021;Zazanashvili et al., 2011), a pattern already observed in Azerbaijan (Buchner et al., 2020). Among the Caucasian broadleaved trees, the beech forests are potentially the highly threatened communities -the current distribution of the species may be reduced by over 45% in this century and largely disappear in Azerbaijan and Armenia (Zazanashvili et al., 2011). The most recent projections are far more pessimistic, indicating only limited suitable areas in the North Caucasus and Iran (Dagtekin et al., 2020;Khalatbari Limaki et al., 2021). Additionally, the lower-elevation populations of this species may be at higher drought risk as was indicated by dendroclimatic analysis (Martin-Benito et al., 2018). This expected reduction would threaten the stability of the forest in the middle mountain belt, where beech dominates, leading to a pronounced biodiversity loss in the region. In this context, the recognition of factors that governed climate-driven range shifts of species is needed to assess the vulnerability to future climate threats (Manel & Holderegger, 2013).
Understanding the patterns of microevolutionary responses of tree species to climate changes remains challenging due to the complex factors involved in the process. These include the uncertainty of future climate scenarios, limitations of species distribution projections, and doubts related to the adaptive potential or the spatio-temporal environmental heterogeneity across species ranges (Alberto et al., 2013;Capblancq, Fitzpatrick, et al., 2020).
Nevertheless, recognition of spatial variation in genetic composition can give insights not only into the impact of the past climate on species' biogeography but also on current population dynamics, particularly the possible genetic consequences of range shifts and, in long term, understanding the effects of climate change (Gavin et al., 2014;Hoffmann et al., 2015). Genetic variation across the range of a species is determined by the interplay of demographic, ecological, and evolutionary processes (Manel & Holderegger, 2013;Orsini et al., 2013). Several theoretical patterns can affect the distribution of genetic diversity, including isolation by distance (IBD), isolation by environment (IBE), isolation by resistance (IBR), or isolation by colonization (IBC; Orsini et al., 2013). However, disentangling the relative contribution of geographic, historical, and contemporary landscape factors affecting these patterns is challenging because they are often spatially correlated, leading to overlapping effects (Nadeau et al., 2016;Orsini et al., 2013). Unraveling the processes underlying the spatial genetic patterns and quantifying the importance of environmental variables in structuring population genetic variation is crucial to managing species and ensuring their long-term sustainability in a changing environment (Hoffmann et al., 2015;Manel & Holderegger, 2013;Orsini et al., 2013). This is especially F I G U R E 1 Distribution range of Oriental beech and the major regions of the Caucasus ecoregion.
important for species occurring in heterogeneous mountainous environments, which are particularly sensitive to the impacts of climate change (Beniston, 2003).
To test for factors shaping spatial variation in neutral genetic composition in a highly heterogeneous landscape, we focused on Oriental beech (Fagus orientalis Lipsky), the most ecologically and economically important tree in the Caucasus (Tarkhnishvili, 2014;Zazanashvili et al., 2011). Its current distribution (Figure 1) includes the Northern Anatolian Mts., the Caucasus Mts., the Talysh Mts.
(Turkey), the Colchis (western Caucasus), and the Hyrcanian (Iran) regions were the main refugia for forest trees, including Oriental beech (Connor & Kvavadze, 2009;Dagtekin et al., 2020;Leroy & Arpe, 2007; E. Ramezani et al., personal communication;Shatilova et al., 2011). While the postglacial migration of the Caucasian temperate forest mostly relied on the Colchis refugium (Connor & Kvavadze, 2009;Tarkhnishvili et al., 2012), the Hyrcanian area is perceived more as the sanctuary of the Neogene flora with limited input into the re-colonization (Akhani et al., 2010). Growing evidence highlights the asymmetrical contribution of the Colchis and Hyrcanian refugia in shaping the modern patterns of genetic structure, suggesting west-east postglacial expansions in the Caucasus and the predominant role of the Colchis (Dering et al., 2021;Parvizi et al., 2019;Tarkhnishvili, 2014). The other detected pattern of interspecific divergence in the Caucasus reflects the vicariance process in these isolated glacial refugia (Christe et al., 2014;Maharramova et al., 2018). Therefore, the Caucasus ecoregion offers an excellent abiotic template to investigate the effects of multiple landscape factors on the contemporary genetic structure of Oriental beech.
We focus on conceptual frameworks that point out the interplay of the neutral and adaptive processes in structuring the neutral genetic diversity in species, as proposed by Orsini et al. (2013). However, due to methodological constraints related to using neutral markers, we mainly discuss neutral processes with some indirect hint about adaptive divergence. Based on available studies (Connor, 2006;Dagtekin et al., 2020;Dering et al., 2021;Tarkhnishvili, 2014), we assumed that the current genetic patterns in Oriental beech have mostly been governed by the colonization history but modified by environmental and adaptive processes. Therefore, we expected decreasing genetic diversity away from the main refugial area due to repeated founding events along the migration routes (Hampe & Petit, 2005). On the contrary, the complex landscape of the Caucasus could induce adaptation to specific habitats promoting intraspecific divergence resulting in a detection of the IBE pattern (Orsini et al., 2013). The adaptive processes may also interact with the neutral ones resulting in the IBC pattern when the local adaptation reinforces the founder effects during range expansions and drive IBC under a monopolization scenario. In this case, the founder effect leads to considerable genetic differentiation among populations and no clear link between the genetic pattern and the spatial and environmental gradients (Orsini et al., 2013). In addition, the long persistence of the species in the isolated refugia could drive a Colchic-Hyrcanian genetic split among populations in those subregions. However, given the species' high potential for gene flow, we may expect a partial eroding of the genetic signal left by the historical factors, leading to overall moderate differentiation. Assessing the future persistence of Oriental beech populations in the Caucasus requires understanding which extrinsic factors determined the current patterns of genetic diversity and connectivity while accounting for their complex evolutionary history.
Here, we applied landscape genetics and ecological niche modeling, aiming at disentangling the historical and contemporary processes underlying the neutral genetic structure of Oriental beech across the South Caucasus. Specifically, we address the following questions: (1) Is genetic diversity spatially structured across the landscape, (2) If yes, what historical, environmental, or spatial processes drive detected patterns of genetic diversity and differentiation? and (3) What are the potential risks involved due to possible changes in the species distribution under future climate projections?
By understanding how the species' genetic structure is associated with current climate variables, we can make the first approximations about the potential risks involved under a future climate (Manel & Holderegger, 2013). Finally, we discuss the implications of the results for the conservation and management of Oriental beech, a keystone tree species of forest ecosystems in the Caucasus.

| Population sampling and genotype acquisition
Sampling covered 32 natural populations of Oriental beech (857 individuals) collected over the entire species range in the South Caucasus ( Figure 2; see Appendix S1: Table S1.1; Appendix S2: Figure S2.1).
Specifically, 19 populations were sampled in the Greater (GC) and Lesser (LC) Caucasus, 11 populations in the Azerbaijan part of the Eastern Greater Caucasus (AZ), and two populations in the Talysh Mts. (southeastern Azerbaijan), which represents the Hyrcanian forests (HZ).
The confidence interval for F ST was determined using bootstrap resampling over loci method with 10,000 replications. To measure the extent of differentiation within the regions, and among the populations, the pairwise F ST following Weir and Cockerham was calculated with ENA correction.
The M-ratio method (Garza & Williamson, 2001) implemented in INEST was used to detect the signature of a recent bottleneck. The M-ratio (MR) was estimated by simulating analysis with 10 5 coalescent replicates under the two-phase mutation model (TPM) assuming a proportion of one-step mutations (ps) of 0.22 and a mean size of multi-step mutations (Δg) of 3.1. The significance of the deficiency in the M-ratio was tested using the Wilcoxon signed-rank test.

| Range-wide population structure
To define the population's genetic structure and admixture, we used STRUCTURE v.2.3 (Pritchard et al., 2000). We assumed admixture models and correlated allele frequencies without prior information on population memberships. Ten replicate runs of independent F I G U R E 2 Locations of the sampled populations of Oriental beech in Georgia (GC, Greater Caucasus; LC, Lesser Caucasus) and Azerbaijan (AZ, East Caucasus; HZ, Hyrcania). Spatial distribution of genetic diversity across the landscape (left panel) with the relationship between both genetic parameters and distance from putative LGM refugial area (Dist LGM ; right panel) -a represents allelic richness (A R ) and b expected heterozygosity (H E ). The population abbreviations as in Table 1; Appendix S1: Table S1.1. subsampling were performed for each genetic cluster (K), ranging from one to 33 with a burn-in period of 10 5 steps followed by 2 × 10 5 MCMC iterations. Following Cullingham et al. (2020), we applied different K-selection methods, including the log probability of the data (Ln Pr (X|K); Pritchard et al., 2000), Evannos' ΔK (Evanno et al., 2005) and the algorithm based on the mean or median membership coefficient (Q) (Puechmaille, 2016). To obtain the K-selection plots, we used StructureSelector (Li & Liu, 2018) while CLUMPAK (Kopelman et al., 2015) was used to summarize and visualize the replicate runs.
According to Puechmaille (2016), we considered clustering results in which a mean membership coefficient (Q) given to genetic clusters is >0.5 to exclude the spurious cluster, which constitutes a doubtful biological grouping. However, being aware of the complexity of the K-selection procedures, we included all clustering results that warrant biogeographic interpretation (Cullingham et al., 2020).

| Species distribution modeling
We used species distribution models (SDMs) to calculate the three landscape metrics: current and past climatic suitability, and distance from the climatically stable area in the LGM. As our point was also to predict the possible changes in the species distribution in future, we constructed the theoretical distribution of the species under different future climate scenarios.
The species occurrence acquisition is detailed in Appendix S2, and for the final modeling procedure, the occurrences dataset hosted 810 unique records (see Appendix S2: Figure S2.1).
The maximum entropy approach implemented in MaxEnt 3.4.1 (Phillips et al., 2004) was applied to build the models. To construct the model of the species' potential distribution for current condition and for future projections (2061-2080), a set of 19 bioclimatic variables at 30 arc-sec resolution were retrieved from CHELSA 1.2 (Karger et al., 2017). Further, these data were upscaled to match the resolution and extent of the bioclimatic variables in QGIS. We used the variance inflation factor (VIF) to eliminate predictor collinearity using the vif function implemented in the usdm R package (Naimi et al., 2014). Variables with large VIF values (>5) were excluded one by one using a stepwise procedure. Finally, the resulting dataset contained nine environmental variables: the annual mean temperature (bio1), isothermality (bio3), temperature seasonality (bio4), mean temperature of the wettest quarter (bio8), mean temperature of the driest quarter (bio9), precipitation seasonality (bio15), precipitation of the warmest quarter (bio18), and the precipitation of the coldest quarter (bio19). The same set of bioclimatic variables at 2.5 arcmin resolution obtained from PaleoClim (Fordham et al., 2017) was applied for past projection during the Last Glacial Maximum (LGM; ca. 21 ka). For past projection, data were obtained from PaleoClim (Fordham et al., 2017); for current and future conditions from CHELSA 1.2 (Karger et al., 2017). The distribution of the species during the LGM (ca. 21 ka) was projected using the Community Climate System Model (CCSM4; Karger et al., 2017), while future projections (2050-2080) based on the Coupled Model Intercomparison Project Phase 5 (CMIP5) following the Representative Concentration Pathways -RCP 4.5 and RCP 8.5 scenarios (Collins et al., 2013).
Maxent was run with 100 replicates using bootstrap resampling, the maximum number of iterations was set at 10 4 , and the convergence threshold was set at 10 −5 with the logistic output of the model prediction for suitability. The "random seed" option was applied to validate the models, where 20% of the occurrence points were random sampling as test data, the remaining points were used as training data, and a random test partition was used for each run. Model accuracy was evaluated using the area under the curve (AUC) values of the receiving operator curve (ROC) as a thresholdindependent evaluation metric (Mas et al., 2013). Results of SDMs across the landscape were visualized using QUANTUM GIS 3.24.0 "Tisler" (QGIS.org, 2022), while habitat suitability and average altitude in the theoretical range of the species were calculated in SAGA GIS (Conrad et al., 2015). To illustrate the extent of the future environmental change among major distributional domains of the species (Greater Caucasus, Lesser Caucasus, Azerbaijan, and Hyrcania), we compared the bioclimatic parameters that had the highest contribution to the SDMs and the current and future habitat suitability using bar charts.
To define populations with high priority for conservation based on genetic data, we applied the Reserve Selection algorithm imple- Moreover, to identify the climatic refugia for Oriental beech, we estimated areas of stability (habitat suitability >60%) defined as a region of overlap between the projected future (RCP4.5 and RCP8.5) and current distribution patterns that support the long-term species occurrence in these regions. For this purpose, the binary Maxent model outputs (30 arc-sec) for the future projections were aggregated to the potential current distribution (30 arc-sec) using raster calculator in QGIS. Results of the areas of stability were visualized using QGIS.

| Predictors of genetic diversity and gene flow
To detect the drivers governing the spatial distribution of genetic diversity, we employed generalized linear models (GLMs) to test the hypothesis that past climate may explain the observed pattern (Hampe & Petit, 2005;Hewitt, 2000) using the glm R function (R Core Team, 2022). The hypothesis emphasizes that higher genetic diversity is related to the proximity of populations to LGM refugia and subsequent decrease due to postglacial migration. In the models, we considered five explanatory variables, including current habitat suitability (HS CURR ), distance from LGM refugium in Colchis (Dist LGM ), genetic admixture (G ADMIX ), latitude and longitude, while A R and H E were used as response variables.
To generate the LGM niche centroid of the species distribution, we used habitat suitability predicted by Maxent, applying the "centroids" option in QGIS. The Euclidean distance was used TA B L E 1 Summary of genetic diversity parameters of Oriental beech (AZ, Azerbaijan; GC, Greater Caucasus; HZ, Hyrcanian stands; LC, Lesser Caucasus) and results of the M-ratio test under the two-phase model (TPM) estimated for sampled populations. as a metric of the population distance from the niche centroid.
The extent of admixture based on the STRUCTURE result at K = 2 (see section 3) was estimated using the "genetic admixture index" obtained according to the procedure described by Ortego et al. (2015). Models were compared using the Nagelkerke pseudo Rsquared, Akaike Information Criterion (AIC), and Akaike weights (w i ) calculated using the function compareGLM in rcompanion R package (Mangiafico, 2022) and model.sel in the MuMIn R library (Bartoń, 2020).
A series of distance-based redundancy analyses (dbRDA) were performed to unravel the relative contribution of climate (clim.), geography (geo.), recent migration (mig.), demographic history (anc.), and topographic heterogeneity (top.) in explaining the detected genetic differentiation (Legendre & Legendre, 2012). To do this, correction as the response variable and a set of explanatory/conditioning variables described below. The analyses were run using the function capscale in vegan R package (Oksanen et al., 2020). First, we tested all the different combinations of all explanatory and conditioning variables in the partial dbRDA to define the "pure" effect of variables and ignore the insignificant variables. Specifically, this constrained ordination approach allowed us to decompose the portion of genetic variance explained by each set of the variables and detect the relative effect of a specific variable by removing the confounding effect of the remaining associated variables that can be spatially correlated (Legendre & Legendre, 2012). Significance of the associations was tested using the anova.cca function with 9999 permutations.
To explore the association of genetic composition with local climate within the isolation by environment (IBE) model, we first selected potentially relevant climatic variables to avoid overfitting and collinearity in the subsequent dbRDAs. To do this, we applied the forward selection procedure using the ordiR2step function in the vegan R package based on the significance test with 9999 permutations and the adjusted R-squared. The tested dataset included climatic variables that were identified as exercising a selective pressure on the beech distribution and genetic variation (Capblancq, Morin, et al., 2020;Pluess et al., 2016). These variables include temperature (the maximum "bio5" and the minimum "bio6" temperature), precipitation (annual precipitation "bio12", precipitation of wettest month "bio13" and driest month "bio14"), evapotranspiration (annual potential evapotranspiration "annualPET"), and drought (aridity index "aridityIndex" and relative wetness to aridity "climaticMois- To test for the effect of isolation by distance (IBD), we calculated a matrix of the Euclidean geographical distances among sampled populations (geo.) estimated from a raster layer depicting a "flat" landscape using QGIS. The current migration matrix was estimated using the divMigrate method (Sundqvist et al., 2016) in the diveRsity R package (Keenan et al., 2013).
The circuit theory within the isolation-by-resistance (IBR) model was applied to test the effect of topographic heterogeneity

| Diversity and differentiation
In total, 293 alleles were detected with an average of 22.54 alleles per locus (Appendix S3:

| Population structure
According to ΔK, the most supported number of genetic clusters Although the second-best group was K = 4 ( Figure 3b), which was also supported by the Ln Pr (X|K) method (Appendix S3: Figure S3.3); this was not considered to be as informative. Two of the inferred clusters did not reach a mean threshold value of Q > 0.5 in all populations assigned to given groups, pointing to the presence of a spurious cluster. Consequently, we considered a clustering result of K = 3, which seems to be justified biologically and showed a clear geographical coherence (Figure 3d,e). This revealed a further substructure of the Azerbaijan populations, splitting populations into two groups, roughly consistent with the north-south pattern of differentiation. However, relatively high admixture across the study sites was observed, mostly among the West-Central and Azerbaijan populations, except for the most distinct ones (Figure 3e,f).
All the Hyrcanian stands and also most of the populations from the Greater Caucasus (except for GC_02, CG_08 and CG_11) and Azerbaijan (except for AZ_2 and AZ_07) showed signs of a significant bottleneck (Table 1). However, only three populations from the Lesser Caucasus (LC_06-CL_08) experienced demographic fluctuations.

| Drivers of genetic differentiation
Among all tested GLMs ( respectively), pointing to west-east decreasing patterns (Figure 2).
According to partial dbRDA, migration, topographic heterogeneity, and geographic distance had an insignificant contribution to the structuring of the species' genetic variation (p > .05; Table 3). After excluding these factors, the full dbRDA model including climate, geographic distance, and ancestry produced a strong significant association (adjR 2 = .914; ≤.001), explaining 74% of the total variance ( Table 3). The different partial dbRDA identified that 23% of this explained variance was associated with the pure effect of ancestry (17%; ≤0.001), and climatic variation (6%; ≤0.001). After excluding geographic distance, ancestry and climate significantly explained 73% of the total variation. (Table 3). Genetic ancestry (IBC) still explained the highest proportion of genetic variation, accounting for 17% even when controlling for confounding effects of other variables, while climate variation (IBE) explained 7% of the total variance ( Table 3). The dbRDA plot indicated a significant association of genetic variation with the climatic variables among Hyrcanian (HZ_01 and HZ_02) and most of Azerbaijan (AZ) stands (Figure 4), showing their divergence from the remaining populations. We found that the two first axes explained most of the genetic variance among the populations (81% in total). dbRDA1 was mostly correlated with aridityIndex and annualPET, while dbRDA2 with bio5 and bio13 ( Figure 4).

| Ecological niche modeling
SDMs showed high levels of predictive performance with a similar score of AUC, reaching >0.946. The precipitation of the warmest quarter (bio18) and the temperature seasonality (bio4) were defined as the most important variables limiting the distributional patterns of the species with a relatively high contribution of >68% and >16%, respectively, in all tested models ( Table 4).

The distribution model under the current climatic conditions
properly described the present range of Oriental beech ( Figure 5).
During LGM, most of the Caucasus region was climatically unsuitable for the species, and potentially suitable conditions existed in three main areas: the most eastern part of the Pontic Mts. (Turkey) with the Adjara region (suitability >75%), the Colchis area with Abkhazia, and the adjacent part of Russia (suitability >75%), and to some extent the Hyrcanian region (suitability 40%). Apart from those refugial areas, the Iori Plateau in southeastern Georgia and the north-western part of the Greater Caucasus in Azerbaijan seemed to offer suitable habitats for the species with a relatively high suitability score reaching 40-75%. Furthermore, the residual areas in the Ganja-Gazakh and Karabakh regions in Azerbaijan were also indicated as climatically suitable for the species but with lower support F I G U R E 3 Spatial genetic structure estimated for the Oriental beech populations across the South Caucasus based on nSSRs using STRUCTURE for K = 2 (above) and K = 3 (below). Pie charts represent the genetic ancestry of each population across the study site (a and d). Admixture assignment of each individual to the inferred K clusters was visualized as barplots; each bar denotes the individual proportion of each of the detected genetic lineages (c and f). K-selection plots according to Evanno's method (2005, b) and Puechmaille (2016, e) approaches show the highest value at K = 2 and K = 3 as the most likely number of clusters, respectively. Population abbreviations as in Table 1; Appendix S1: Regarding the future predictions under RCP4.5, significant changes in habitat suitability in the Hyrcania and western part of the species range are not expected ( Figure 5). Nevertheless, a reduction

F I G U R E 6
Populations of Oriental beech with priority for conservation inferred with reserve selection algorithm implemented in DIVA-GIS in relation to the potential future distribution under the pessimistic scenario (RCP8.5, ca. 2070). Disk diameters are proportional to the value of genetic parameters, following the figure legends (left panel). Bar charts presenting bioclimatic variables with the highest contribution in SDMs, variables significant association with genetic structure, habitat suitability and average shifts in elevation for the current and future projections (right panel). Population abbreviations as in Table 1; Appendix S1: Table S1.1. (Christe et al., 2014;Ibrahimov et al., 2010;Maharramova et al., 2018;Tarkhnishvili, 2014).
Given the detected significant eastwards gradient of genetic diversity, we can conclude that the Colchis refugium likely acted as the major source of postglacial expansion. Notably, the west-east genetic pattern seems to be also consistent with the assumptions of the leading-edge model of range expansions related to densitydependent processes (Hampe & Petit, 2005). This pattern has been previously documented for trees that have undergone a longdistance migration pattern (Roberts & Hamann, 2015). Given the general assumptions regarding the refugial area (Hewitt, 2000), the high allelic diversity and private alleles found in populations from the long-distance pollen dispersal, distances even up to 1000 km, is frequent in beech (Belmonte et al., 2008;Piotti et al., 2012). A similar pattern of the efficient pollen-mediated gene flow among distinct refugia during postglacial expansion has been shown for Abies alba (Liepelt et al., 2002;Piotti et al., 2017) and Pinus banksiana (Godbout et al., 2010). Wider sampling in the East Caucasus would shed light on the postglacial migration in this area.

| Drivers of genetic and differentiation patterns
Considering the postglacial history of the Oriental beech mostly related to single Colchis refugium, and environmental gradients present in the study area, it seems that our inference is burdened by uncertainty due to the correlation of spatial-environmental factors.
However, we were able to disentangle the forces structuring neutral genetic diversity across species' ranges, applying the variance partitioning approach that reduces the confounding effects of potential spatially correlated predictors and quantifies its relative influence (Legendre & Legendre, 2012).
Considering the topographic complexity of the Caucasus and the IBC hypothesis that involves both historical and adaptation processes as drivers of the population's divergence (Orsini et al., 2013), we expected to find strong geographically structured diversity with a clear split between the Greater and Lesser Caucasus. However, we found no support for topographic complexity being an important factor in genetic structure. Genetic distinctiveness between these mountain ranges for another wind-pollinated tree, Pinus sylvestris (Dering et al., 2021), has been explained by the direction of prevail-

The clear west-east gradient of genetic differentiation in
Oriental beech could suggest a strong pattern of IBD. However, given that IBD does not account for the landscape heterogeneity (Jenkins et al., 2010) and can interfere with the alternative patterns of population structure resulting from colonization history and landscape resistance to gene flow (Orsini et al., 2013;van Strien et al., 2015), this seems to be an unrealistic scenario due to the oversimplification of processes involved. Accordingly, after controlling for the confounding effect of genetic ancestry and climate, the geographic distance by itself had comparatively little contribution to the observed genetic pattern (1% of the total variation, p > .05). It seems that the mid-elevation areas in the west-central In the absence of palaeobotanical evidence for a cryptic refugium in the eastern part of the South Caucasus that could act as a source of eastward colonization, we cannot conclusively state whether the distinct genetic composition of Azerbaijan populations is consistent with the IBC or IBE. Since postglacial colonization is able to generate patterns similar to IBE (Hampe & Petit, 2005;Orsini et al., 2013), it seems that detected divergence is a result of recent postglacial history rather than vicariance process in an isolated cryptic refugium. Additionally, we did not find the accumulation of private alleles and high genetic variation in the Azerbaijan populations, which could support the presence of a cryptic refugium in this area. Conversely, the initial pattern produced by vicariance might have been partly swamped by the current relatively high gene flow among west-central Caucasus and Azerbaijan populations resulting in low genetic differentiation (F ST = 0.011). Consequently, the current distinctiveness of the Azerbaijan sites could be a weak signal of the initial founder effect originating from the colonization stage. According to IBE model, the selection against maladapted migrants may allow the genetic signal of the initial structure to be preserved in the neutral diversity for generations (Orsini et al., 2013).
The divergence of the Hyrcanian and most of the Azerbaijan stands (AZ_04-AZ_11) could also have been caused by local adaptation, given the climatic distinctiveness of the East Caucasus. The detected significant signal of IBE, accounting for 6% of the total variation, suggests that the genetic composition is partially structured by local climate. Specifically, aridity, maximum temperature, precipitation of the wettest month, and annual potential evapotranspiration were significantly associated with genetic distance. According to the autecology diagram, these range-edge populations can be considered as ecologically marginal ( Figure 6). Such a distributional pattern implies the development of local adaptations. However, due to methodological constraints, our results are not a pertinent proxy of adaptive divergence, which requires the detection of genomic signals of adaptation. Nevertheless, the selectively neutral markers may show some association with the environment due to genome hitchhiking leading to the IBE patterns (Nosil et al., 2008), which means that the hypothesis on the contribution of local adaptation to the genetic structure of Oriental beech remains valid. Indeed, a strong association between neutral genetic composition and environmental gradients has been found in other tree species (Muniz et al., 2022;Sork et al., 2010).

| Conservation implication
In contrast to some studies (Dagtekin et al., 2020;Khalatbari Limaki et al., 2021), our SDM models are not so pessimistic about the future theoretical distribution of Oriental beech, especially in Turkish and Hyrcanian parts of the range. Nevertheless, much of the currently highly suitable areas for the species may be lost. The most prominent changes are the distributional contractions projected in the Azerbaijan part of the Greater Caucasus, Armenia, and eastern Georgia ( Figure 5). Moreover, the range shifts westward and may show a twofold increase in elevation under the most pessimistic scenario (Figure 6). The shifting to the higher elevation of the species can be mostly explained by temperature increases because other climatic trends (e.g., precipitation) are not generally related to elevation. According to the climate projection, the mean temperature in the Caucasus Mts. is expected to rise by at least 3°C by the end of this century compared with the current condition. Higher temperatures are assumed to increase the intensity of soil drought due to the forcing effect on potential evapotranspiration (Bergh et al., 2003).  (Gavin et al., 2014;Karger et al., 2017). To reduce this uncertainty, we used climatic data from CHELSA that has higher accuracy in mountain-specific conditions (Brown et al., 2018;Karger et al., 2017). Additionally, using the occurrence dataset drawn from map grid cells, as was done in Dagtekin et al. (2020), can be a source of model bias (Konowalik & Nosol, 2021).  (Koskela et al., 2013). They can be directly applied also in the Caucasian populations of Oriental beech. Our results provide additional information regarding genetic diversity that can support the process of GCUs establishment. Moreover, in light of the assumption that species adaptation to climate change mostly relies on the standing genetic variation , special attention should be paid to the population from the Adjara-Imereti Range. These populations host the highest and most unique neutral genetic diversity that is of crucial conservation and management priority and can be highly relevant for the future resilience of the species. Additionally, populations from the Trialeti Range (LC_06-LC_08), Khevsureti (GC_06), Gombori Range (GC_08), and Azerbaijan stands (AZ_01 and AZ_05) should be preserved in the context of maintaining a high spectrum of genetic diversity needed for sustainable beech forest management.
However, given the adaptive potential of range-edge populations to climate change (Fady et al., 2016;Hampe & Petit, 2005;Rehm et al., 2015) and that the species ability to persist under such changes will be determined by the responses of the local populations (Aitken et al., 2008), the Oriental beech populations at the range-edge should also be considered. Our SDM showed that the eastern Caucasian gene pool of the species is expected to be seriously vulnerable because of increases in temperature and aridity ( Figure 6), especially the peripheral Azerbaijan populations that already occur in marginal conditions and display low gene diversity.
The detected excess of inbreeding and signs of bottlenecks may suggest that adverse demo-genetic processes are already present in these populations. On the other hand, these populations may potentially harbor important adaptive properties generated under such environmental constraints (Aitken et al., 2008;Rehm et al., 2015). Given the projected extreme decreases in precipitation in the eastern domain of the species range, the probable intensification of the stochastic genetic processes may pose a risk to that unique gene pool. Another possible consequence of the climate-driven range shifts might be the loss of landscape connectivity, triggering strong genetic drift. Furthermore, the detected strongly asymmetric gene flow among the Georgian and Azerbaijan population may also have serious evolutionary consequences related to adaptation lags of range-edge populations due to receiving maladaptive alleles (Aitken et al., 2008;Fréjaville et al., 2020).

| CON CLUS I ON S AND FUTURE DIREC TIONS
We are aware that a complete understanding of how ecologically marginal populations of Oriental beech may cope with climate change adaptation requires a detailed investigation including a genome-environmental association approach to identify a signature of local adaptation and the recognition effect of gene flow on adaptation (Capblancq, Fitzpatrick, et al., 2020). The SDM assumes a genetic homogeneity, and incorporating the adaptive genetic variation in climate change vulnerability assessment could deliver more reliable projections. Our study is the first approximation of the potential risks involved in climate change and induces far-reaching thinking about the need of applying management solutions dedicated to maintaining the genetic resources of Oriental beech (Fady et al., 2016).
This study enriches our understanding of the evolutionary history of Oriental beech and the forces that shape its neutral genetic composition in the South Caucasus. Nevertheless, several other questions remain unanswered, waiting for comprehensive sampling across the whole species range and implementation of more relevant landscape genomic and demographic approaches. We would like to know what is the adaptive genetic potential of Oriental beech, how is it distributed across the species' range, and how this can be helpful for the species in tracking future climate change. These issues are important because of the potential range reduction of one of the most valuable Caucasian tree species, with implications for forest management in Europe (Brang et al., 2016).

ACK N OWLED G M ENTS
The study was supported by the National Science Centre (Grant Number 2017/26/E/NZ8/01049), the Institute of Dendrology, Polish Academy of Sciences and the Poznań University of Life Sciences.
We thank Berika Beridze, Georgian and Azerbaijan foresters for assistance in sampling, Łukasz Dylewski for help in the statistical analyses, and Thibaut Capblancq for his comments. We confirm that experimental design including sampling of the plant material used in this work complied with institutional, national, and international guidelines and legislation.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study, including genotypes, raw dataset used for landscape genetics analyses and